Reduced basis catalogs for gravitational wave templates 
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We introduce a reduced basis approach as a new paradigm for modehng, representing and search- 
ing for gravitational waves. We construct waveform catalogs for nonspinning compact binary co- 
alescences, and we find that for accuracies of 99% and 99.999% the method generates a factor of 
about 10 — 10* fewer templates than standard placement methods. The continuum of gravitational 
waves can be represented by a finite and comparatively compact basis. The method is robust under 
variations in the noise of detectors, implying that only a single catalog needs to be generated. 



Introduction. The second generation of earth-based 
gravitational wave detectors, such as Advanced LIGO 
and Virgo, will become operational in 2014-2015. These 
detectors are expected to directly measure gravitational 
waves (GWs), with likely event rates of 0.4 — 400 per year 
for binary neutron stars (NS) and 0.4 — 1, 000 for binary 
black holes (BH) [1] . Direct detections would allow tests 
of general relativity in the nonlinear regime as well as ac- 
cess to portions of the universe otherwise unobservable. 

Compact binary coalescences (CBCs), which consist of 
a pair of NSs and/or BHs inspiraling and merging, are 
considered to be one of the most promising sources of 
gravitational waves. The preferred method to search for 
GWs from CBCs is to employ matched filtering, which 
compares data from a detector to a bank of possible tem- 
plate waveforms and checks for a strong correlation be- 
tween them. In low mass searches, the inspiral regime 
dominates the observable signal but as the mass increases 
the merger regime becomes increasingly relevant. The 
merger regime requires numerical simulations ~ even if 
used only for calibration of semianalytic models. Given 
the number and cost of these simulations, knowing an 
optimal choice of parameters is critical in order to limit 
the number of large simulations needed to accurately rep- 
resent the variation over the parameter space. For this 
reason, it is desirable to seek a method that builds a tem- 
plate bank of waveforms by sequentially selecting only 
the most relevant points in the parameter space. 

In addition, once a waveform catalog is constructed, 
there is a significant computational cost in performing 
an actual search for GWs due to the size of the cata- 
logs. Real-time analysis of the data is critical to gener- 
ate alerts to search for electromagnetic counterparts and 
enable multimessenger astronomy [2 -4J. With the stan- 
dard catalog placement method the number of templates 
needed grows rapidly with the dimension P of the pa- 
rameter space (as (1 — MM)^-^/^, with MM the minimal 
match [5]) and such an approach could become imprac- 
tical for searches of spinning binaries or other complex 



physical systems. 

Reduced Basis Method. The RB framework [6| con- 
structs a global basis rather than using local methods and 
can be seen as an application-specific spectral expansion. 
In such an approach one seeks to enable a rapid online 
evaluation of the reduced model at the expense of having 
to build the basis prior to the application. While such 
an approach is not suitable for certain types of applica- 
tions, it is particularly well suited for those requiring near 
real-time or many inquiry response as is the case in the 
present application. It has the following advantageous 
features over standard model reduction techniques such 
as Proper Orthogonal Decompositions, Singular Value 
Decompositions, or Principal Component Analysis (see 
J\ for a general review of these methods and [51 12] for 
applications to GWs), see also [TU] : 

1. It is applicable to situations in which one must 
choose the most relevant parameters on the fly. 

2. It yields nested, hierarchically constructed cata- 
logs which can be easily extended. If Cn — 
{hi, ... ,hiy} is a catalog from the RB method then 
adding additional waveforms for higher accuracy 
implies that the resulting catalogs contain the pre- 
vious ones, Cn C Cjv+i C Cn+2 ■ ■ ■ ■ 

3. It is highly computationally efficient. The cost of 
adding a new member to an existing catalog of size 
N is independent of N. Hence, the total cost of 
generating a catalog of size N scales linearly with 
N, in contrast to many other approaches. 

4. It yields catalogs that are nearly optimal in terms 
of the error in approximating the whole spectrum 
of GWs by a compact set of basis elements. This 
error is measured in the L°° norm, ensuring a strict 
upper bound over the entire parameter space. 

A gravitational wave is a function of time (or fre- 
quency in Fourier space) and of the P parameters /Z = 
{^1, . . . , fip} associated with the source. We denote each 
of them simply by hp and do not explicitly write the time 



2 



or frequency dependence. Let % be the space of all nor- 
malized GWs for the considered source(s). Although % 
is a not a linear space (the sum of two waveforms is not 
a waveform), we show that it can be represented by a 
linear one with arbitrarily high accuracy. We start with 
a theoretical description of our approach, followed by a 
description of an actual implementation. 

We are interested in approximating % by the best lin- 
ear combinations of members = hp^p. of a catalog 
Cn = {'^i}iLi- AH such linear combinations form the 
reduced basis space Wn = span (Cat). The waveforms 
that make up this catalog could be optimally chosen so 
that the error in representing Ji with Wn is minimized 
over the choice of TV catalog members. Such an optimal 
error is given by the Kolmogorov A^- width [TT] , 

(iAr('H) = min max min (1) 

Cn P u£Wn 

That is, one computes the error in the best approxima- 
tion of hp by a member of Wn, then finds the parameter 
fl yielding the largest error, and lastly finds the small- 
est such error for all possible A^-member catalogs. Here, 
the norm in Eq. ([l]) is calculated from the complex inner 
product (•, •), which is related to the standard overlap of 
Wiener filtering by 43?[(-, •)], such that for two waveforms 
F and G in Fourier space. 
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where <S'„(/) is the power spectral density (PSD) of the 
detector. 

Finding a catalog that exactly achieves the A^-width 
is a computationally demanding optimization problem. 
Instead, we use a greedy approach, which is an inexpen- 
sive and practical procedure for hierarchically generating 
catalogs that nearly satisfy the A^- width [12]. 

One constructs a catalog by first choosing a waveform 
for an arbitrary parameter value. A basis vector ei is 
then identified with this waveform, ei — hp-^ , and the 
catalog is Ci = {vPi = hp^ }. To add another waveform to 
the catalog, one seeks the parameter value jl2 that max- 
imizes \\hp — Pi{hp)\\ where Pi{hp) = ei{ei,hp) is the 
(orthogonal) projection of hp onto Wi. We call this step a 
greedy sweep. The waveform corresponding to /22 is added 
to the catalog so that C2 = {^'i, ^'2}. The new basis vec- 
tor 62 is then constructed via Gram-Schmidt orthonor- 
malization. Notice that Ci C C2, which demonstrates 
the hierarchical nature of the catalogs generated. Addi- 
tional members of the reduced basis catalog are generated 
by mathematical induction. At each step one picks up the 
parameter value jlj that maximizes \ \h{^) — Pj^i{h{fl))\\. 

It can be shown [T^] that if the decay of the A^-width 
with A^ can be bounded by an exponential, 
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FIG. 1: Error in approximating the space of waveforms by 
a discrete catalog for BNS inspirals with Initial LIGO. For 
reduced basis, the error is the square of the greedy error Q 
while for metric placement the error is (1 — MM) with MM 
the minimal match. The lower panel shows the extrapolation 
of the maximum number of RBs generated for an infinitely 
large training space. The fit shown (red) excludes the two 
points with largest x, which change the asymptotic value by 
0.2. 



for some real c and a, then the decay of the maximum 
error for a catalog Cn generated by this approach, which 
we call the greedy error En, is also exponential. 



e^v = max 
P 



PN{hp) 



(3) 



where PN{hp) = '^f^ii^i, hp)ei and d, /3 depend on c, a 
(see [12] for more details). Similar results hold in the 
case of power-law fall-off, i.e. dN{T~L) < BN^'^ implies 
En < BN~''. Note that ejv is a bound on the error 
between a waveform and its representation, and that 
= maxp{l — 3i[{hp,PN (hp))]), so that e]^ is an er- 
ror on the overlap directly comparable to (1 — MM). 
Given that GWs appear to depend smoothly on the pa- 
rameters jl, we expect dNCH), and hence the greedy 
error en, to decay rapidly (in fact exponentially) with 
N, which is a key feature of this method. Notice that 
([3]) implies that any waveform can be represented as 
hp = PN{hp) + 5hp{f) where \\5hp{f)\\ < En- There- 
fore, if En is of the order of numerical round-off then, in 
practice, the projection of hp onto Wn equals the wave- 
form itself. In addition, the number of RBs needed to 
represent any hp is comparatively small (see below). 

In the case in which one is numerically solving equa- 
tions on the fly while building the RB, the error ([s]) is re- 
placed by an inexpensive error estimate evaluation (such 
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as a residual), which is referred to as the weak greedy 
approach. Once the next parameter value is chosen one 
solves the full problem for it (referred to as the offline 
stage) and proceeds to the next greedy sweep. In any 
greedy approach, the maximum over fl is searched for, 
in practice, using a training space of samples fl. How- 
ever, since this is done as part of the offline process, the 
training space can be finely sampled and one can take ad- 
vantage of the observation that evaluations for different 
parameters values are decoupled and, hence, embarrass- 
ingly parallel. 

If one attempted a matched filter search with a RB 
catalog Cat by filtering each basis function against the 
data and maximizing over arbitrary linear combinations 
of these filter outputs, one would of course get a very 
high false alarm rate. Instead, it is important to al- 
low only linear combinations that correspond to physical 
waveforms. To be more specific, one could easily store 
the matrix of overlaps between waveforms in the origi- 
nal template bank (the training space) and the reduced 
basis, i.e., aij ~ {ei,hp.). In fact, our algorithm pro- 
vides this reconstruction matrix as output. A matched 
filtering computation may then be performed by inte- 
grating the incoming signal s against each member of 
the basis (s^Ci) and using the reconstruction matrix aij 
to recover the matched filtering integral of the signal 
with any template hp. in the original bank. Explicitly, 
{s,hjx-) = '^i{s,ei)aij. In this way, using the reduced 
bases is equivalent to using the original waveform space 
but with many fewer matched filtering integrals to com- 
pute for a given signal. Hence, using RB yields no in- 
crease in the false alarm rate. 

Catalogs for compact binary inspirals. We discuss our 
results for constructing reduced bases for "chirp" grav- 
itational waveforms for binary inspirals without spins 
[H] . We use the 2nd order post- Newtonian accurate 
waveforms in the stationary phase approximation, which 
are known in closed form, so that the parameter space 
is two-dimensional (the binary's masses). For simplicity, 
we take the coalescence time and phase to be constant 
for each waveform. 

Fig.[l]shows results for the greedy error using a reduced 
basis model for inspirals of binary neutron stars (BNS) 
with mass components in the range [1-3] (for Initial 
LIGO with a lower frequency cutoff at 40 Hz) compared 
with the standard metric template placement method [5] . 
After a slowly decaying region, the reduced basis model 
gives very fast exponential convergence decay, which can 
be fitted by e% = ae-''^" with a = 9.65 x 10"'', b = 0.598, 
p = 1.25. The metric method yields approximately linear 
decay for a two-dimensional parameter space. As already 
mentioned, this decay becomes slower as the dimension- 
ality P of the parameter space increases. The fast decay 
of the reduced basis model allows a representation of the 
whole set of gravitational waves for these sources and 
mass ranges to within machine precision. We have found 
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FIG. 2: The points show the parameter values chosen for the 
catalog of BNS and Initial LIGO. The density of parameter 
values is shown using a coloramp as well as histograms. 

the same feature in all mass ranges that we have explored. 
This leads to the rather remarkable finding that for all 
practical purposes the set of relevant gravitational wave- 
forms in compact parameter regions appears to be finite 
dimensional. When increasing the number of samples x 
in the training set we find the following fit for the number 
of RB for machine precision error, N = a + bx~^^'^ +cx~^ 
with a = 921, b = -2090, c = -9.18 x 10^ for the case of 
Fig. [l] In particular, in the limit x — > oo only 921 bases 
are needed to represent, within numerical accuracy, the 
full space of waveforms H for this range of masses for 
BNS inspirals. 

Figure [2] shows the chosen parameter values in the 
chirp mass vs symmetric mass ratio plane and a density 
plot of the number of RBs. The histograms highlight that 
most values are picked for (nearly) equal mass systems 
of low chirp mass. 

Table |T] shows the number of RB that we need to rep- 
resent, for different overlap error tolerances, inspirals of 
BNS and stellar size binary black holes (BBH, with mass 
components in the range [3-30] A/q). The limit x — > oo is 
not taken here for simplicity so the RB values listed in 
Table |T] are slightly underestimated. 

Sensitivity to Nonstationary Noise. The PSD of any 
ground-based interferometer will fluctuate in time due to 
changes in environmental noises and other factors. Since 
the PSD weights the inner products used to construct the 
reduced basis, one might worry that a new RB needs to 
be constructed for any variation in the PSD. 

Remarkably, we find indications that the RB con- 
structed assuming a fiducial PSD is highy robust against 
rather large perturbations. From a histogram of the sen- 
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TABLE I: Number of reduced bases/templates for different 
target accuracies with the reduced basis (RB) and template 
metric (TM) approaches for binary neutron stars (BNS) and 
binary black holes (BBH), using spinless chirp waveforms. We 
assume a lower frequency cutoff of 40 Hz for Initial LIGO and 
10 Hz for Advanced LIGO and Virgo. The error is given by 
e% for RB and (1 - MM) for TM. 



Detector 


Overlap 
Error 


BBH 


BNS 


RB 


TM 


RB 


TM 


InitLIGO 


10"^ 
10-'^ 
2.5 X 10"^^ 


165 
170 
182 


2,450 
1.2 X 10'' 
5.9 X 10^^ 


898 
904 
917 


10,028 

4.3 X 10'' 

1.4 X 10" 


AdvLIGO 


10"^ 
10"^ 
2.5 X 10"" 


1,058 
1,687 
1,700 


19,336 
1.5 X 10^ 
2.3 X 10" 


5,395 
8,958 
8,976 


72, 790 
4.9 X lO'^ 
5.6 X 10" 


Adv Virgo 


2.5 X 10"" 


1,395 
1,690 
1,703 


42, 496 
3.1 X lO'^ 
4.8 X 10" 


7,482 
8,960 
8,977 


156,127 
8.3 X lO'^ 
6.0 X 10" 



sitivity of the LIGO interferometers during a portion of 
LIGO's fifth science reported in [T5], we conclude that a 
20% increase or decrease in the sensitivity to BNS sig- 
nals is a large but realistic fluctuation to the sensitivity 
of the Initial LIGO interferometers. Therefore, we con- 
structed smooth deformations of the Initial LIGO design 
PSD meant to simulate variations in the seismic, thermal 
and shot noise levels which yield roughly a 20% increase 
or decrease in sensitivity to BNS signals. We find that 
the RB generated for the Initial LIGO design PSD with 
a greedy error tolerance of En can represent any inspi- 
ral waveform within the perturbed PSDs considered here 
with an error of no more than I.Sen- This result implies 
that one needs to compute only a single reduced basis for 
a given source for a particular detector. This is unlike 
the current operating procedure in which a new template 
bank is generated every ^ 2048 seconds because of the 
drifts in the nonstationary noise. We will provide further 
details in a forthcoming paper. 

Conclusions. We have considered the development 
and use of a reduced basis method to template bank con- 
struction and found rapid exponential convergence of the 
waveform catalog over the full parameter space. The cat- 
alog is computationally cheap to derive, hierarchical (i.e. 
if a more accurate catalog is required, elements can be 
added) , can be extended for a computational cost that is 
independent of N, and is robust under changes in a detec- 
tor's noise. We have found that the space of gravitational 
waveforms considered in this paper is essentially finite- 



dimensional for any finite range of physical parameters, 
and conjecture that it is in general the case. Elsewhere, 
we will present a more detailed description of these re- 
sults and further applications of the RB framework. 
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